The following R packages are used in the analysis.

library(tidyverse)
library(DataExplorer)
library(brms)
library(rstanarm)
library(ggdist)  
library(skimr)
library(Amelia)  
library(caret)
library(rpart.plot)
library(gridExtra)
library(lme4)
library(BAS)
library(latex2exp)

# Set theme for ggplot
mytheme = theme_bw()
theme_set(mytheme)

Read the data

d <-  read_csv("PerturbedCyclistExperimentalDataTUDelft2022.csv", 
        col_types = cols(participant_ID = col_character(), 
        outcome = col_factor()))

We reorder participant_ID according to age, filter out initial search data and remove some columns that we don’t further need in the analysis.

d <- d %>% 
  mutate(participant_ID = fct_reorder(participant_ID, age)) %>%
  filter(initial_search==FALSE, exclude=="NO") %>% 
  select(-c(pull_force, comments, initial_search, exclude)) 

d %>% glimpse()
## Rows: 928
## Columns: 19
## $ participant_ID    <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ gender            <chr> "female", "female", "female", "female", "female", "f…
## $ age               <dbl> 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, 61, …
## $ weight            <dbl> 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6, 83.6…
## $ length            <dbl> 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61, 1.61…
## $ skill_performance <dbl> 81.30124, 81.30124, 81.30124, 81.30124, 81.30124, 81…
## $ skill_effort      <dbl> 3.798331, 3.798331, 3.798331, 3.798331, 3.798331, 3.…
## $ velocity          <dbl> 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, …
## $ pull_ID           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
## $ pull_direction    <chr> "CCW", "CCW", "CW", "CW", "CCW", "CCW", "CW", "CCW",…
## $ outcome           <fct> 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0…
## $ angular_momentum  <dbl> 6.7600990, 9.2634520, 11.9062300, 2.9413360, 8.27263…
## $ reaction_time     <dbl> 0.0045, 0.0190, 0.0760, 0.0990, 0.0490, 0.0655, 0.03…
## $ phi_zero          <dbl> -1.176674000, -0.978592300, -1.416520000, -0.9767104…
## $ phi_dot_zero      <dbl> 0.194354000, 0.163115900, -0.013537920, 0.080976880,…
## $ delta_zero        <dbl> -0.030918070, -0.017670920, -0.018878270, -0.0093245…
## $ delta_dot_zero    <dbl> 0.202679100, -0.127448200, -0.195207200, 0.156655100…
## $ psi_zero          <dbl> -0.027438300, -0.003981016, 0.023122360, -0.01823309…
## $ Q_zero            <dbl> -46.488590, -39.377390, 34.093770, 31.780850, -1.862…

Exploratory data analysis

Missing data

d %>% skim()
Data summary
Name Piped data
Number of rows 928
Number of columns 19
_______________________
Column type frequency:
character 2
factor 2
numeric 15
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
gender 0 1 4 6 0 2 0
pull_direction 4 1 2 3 0 2 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
participant_ID 0 1 FALSE 24 22: 60, 2: 60, 15: 60, 11: 60
outcome 0 1 FALSE 2 0: 484, 1: 444

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 0 1.00 43.74 17.67 20.00 29.00 34.00 62.00 75.00 ▇▃▁▅▃
weight 0 1.00 77.01 11.70 58.30 67.90 76.40 85.10 103.30 ▆▇▇▅▂
length 0 1.00 1.81 0.09 1.61 1.73 1.83 1.87 1.98 ▃▃▅▇▃
skill_performance 20 0.98 65.45 26.69 22.32 42.91 63.77 85.70 135.02 ▇▇▇▅▁
skill_effort 0 1.00 4.78 3.56 0.00 2.57 4.18 6.42 18.82 ▇▇▂▁▁
velocity 0 1.00 12.12 4.22 6.00 12.00 12.00 18.00 18.00 ▃▁▇▁▅
pull_ID 0 1.00 22.49 14.72 1.00 10.00 20.00 33.00 60.00 ▇▇▆▃▂
angular_momentum 0 1.00 13.95 5.22 0.00 10.37 13.80 17.54 29.23 ▁▆▇▅▁
reaction_time 17 0.98 0.03 0.02 0.00 0.01 0.03 0.05 0.15 ▇▇▂▁▁
phi_zero 24 0.97 -0.98 1.14 -4.41 -1.85 -0.92 -0.13 4.10 ▁▇▇▁▁
phi_dot_zero 15 0.98 0.03 0.14 -0.92 -0.01 0.03 0.07 1.04 ▁▁▇▁▁
delta_zero 24 0.97 0.00 0.03 -0.15 -0.02 0.00 0.02 0.17 ▁▂▇▁▁
delta_dot_zero 15 0.98 0.07 0.41 -2.37 -0.11 0.07 0.25 2.37 ▁▁▇▁▁
psi_zero 24 0.97 0.00 0.02 -0.09 -0.01 0.00 0.01 0.12 ▁▅▇▁▁
Q_zero 15 0.98 13.38 64.76 -320.51 -22.80 11.12 47.29 240.67 ▁▁▇▆▁ 
plot_missing(d, missing_only=TRUE, ggtheme =mytheme)

d %>% missmap(y.labels=d$participant_ID, 
              main="Missingness map, (participant_ID on vertical axis)")

Histograms of predictors on a numerical scale

dcontinuous <- d %>% select(age, angular_momentum, delta_dot_zero, delta_zero, length,
                            phi_dot_zero, phi_zero, psi_zero, Q_zero, reaction_time, 
                            skill_effort, skill_performance, weight)
dcontinuous %>% plot_histogram(ggtheme=mytheme)

Zoom into peak near zero for reaction time

d %>% drop_na() %>%
  mutate(`log10(reaction_time)`=log10(reaction_time)) %>%
  ggplot(aes(x=`log10(reaction_time)`)) +
  geom_histogram(bins=50,colour="white")

Investigate correlations

plot_correlation(dcontinuous, ggtheme = mytheme)

Hence, only length and weight have moderate correlation.

Compute fraction of falling for each velocity and participant

fracfalling <- d %>% mutate(velocity_=as.factor(velocity)) %>%
  group_by(participant_ID, velocity_) %>%
  summarise(n=n(), 
  fraction_falling=mean(outcome=="1"),
  skill_performance=mean(skill_performance),
  participant_ID = first(participant_ID),
  age=mean(age),
  gender=first(gender),
  skill_effort=first(skill_effort),
  average_reaction_time=mean(reaction_time))

fracfalling
## # A tibble: 47 × 9
## # Groups:   participant_ID [24]
##    participant_ID velocity_     n fraction_falling skill_performan…   age gender
##    <fct>          <fct>     <int>            <dbl>            <dbl> <dbl> <chr> 
##  1 9              6            20            0.4               34.3    20 female
##  2 9              12           20            0.5               44.7    20 female
##  3 16             12           20            0.45              60.4    24 female
##  4 16             18           19            0.579             36.0    24 female
##  5 19             12           20            0.45             115.     27 female
##  6 19             18           20            0.45              86.3    27 female
##  7 22             6            20            0.35              46.3    27 male  
##  8 22             12           20            0.5               75.0    27 male  
##  9 22             18           20            0.5               58.0    27 male  
## 10 6              12           20            0.55              85.7    27 male  
## # … with 37 more rows, and 2 more variables: skill_effort <dbl>,
## #   average_reaction_time <dbl>

Visualise

fracfalling %>% 
  ggplot(aes(x=participant_ID, y = fraction_falling, colour=velocity_)) +
  geom_point()

Influence angular momentum

d %>% 
  ggplot(aes(x=angular_momentum, y=outcome)) +
  geom_jitter(height=0.2)

We can make separate panels for each participant and colour according to velocity

d %>% mutate(velocity=as.factor(velocity)) %>%
  ggplot(aes(x=angular_momentum, y = outcome, colour=velocity)) + 
  geom_jitter(height=0.1,size=0.5) + 
  facet_wrap(~participant_ID,nrow=4) + 
  labs(x="angular momentum", y="y")+
  theme(legend.position = "top")

Some QQ-plots

We investigate in which sense the variables phi_zero, phi_dot_zero, delta_zero, delta_dot_zero, psi_zero and Q_zero differ over the binary outcome.

dmeas <- d %>% select(outcome, phi_zero, phi_dot_zero,
                 delta_zero, delta_dot_zero, psi_zero, Q_zero) %>%
                 pivot_longer(cols=contains("zero"), names_to="measurement", values_to="value")
  
dmeas %>%  
  ggplot(aes(sample=value,colour=outcome)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~measurement, scales="free") +
  labs(x="theoretical", y="sample")

We observe similar patterns for all of these 6 variables.

Relation phi_dot_zero and pull direction

d %>% ggplot(aes(sample=phi_dot_zero,colour=outcome, shape=pull_direction)) +
  stat_qq(size=2)

Classification tree and random forest

Classification tree

We fit a classification tree to get a first impression of variables that matter for predicting the outcome variable. The advantage of such a method is that it is insensitive to noise variables that don’t have predictive power.

tr <- train(outcome ~ ., data = na.omit(d), method='rpart', tuneLength=10)
rpart.plot(tr$finalModel,extra=101)

ggplot(varImp(tr))

Random forest

Classification trees are highly variable. Therefore, we fit a random forest. As the fraction of rows containing missing values is relatively small, we fit the RF without rows containing missing values.

d_dropna  <- d %>% drop_na()
tr_rf <- train(outcome ~ ., data = d_dropna, method='rf',
               trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
ggplot(varImp(tr_rf))

We train the random forest again, but now such that participant_ID is not automatically converted to dummy variables. This implies that variable importance is computed for participant_ID by itself, and not for individual participants.

X <- d_dropna %>% dplyr::select(-outcome)
Y <- d_dropna %>% dplyr::select(outcome) %>% deframe()
tr_rf2 <- train(X,Y, method='rf',
            trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
p2 <- ggplot(varImp(tr_rf2))

Refit, but now further excluding participant_ID.

Xsub <- X %>% dplyr::select(-participant_ID)
tr_rf3 <- train(Xsub, Y, method='rf',
            trControl=trainControl("cv"), importance = TRUE, tuneLength=10)
p3 <- ggplot(varImp(tr_rf3))

Put the two figures side by side

grid.arrange(p2,p3, ncol=2)

Not surprisingly, if participant_ID is dropped, then length, weight and age get assigned higher variable importance.

Bayesian Model Averaging

Our analysis has two objectives:

First, we standardise (mean centered, scaled to have unit standard deviation) the numerical variables. This helps to interpret the magnitude of the estimated coefficients.

preProcValues <- preProcess(d, method = c("center", "scale"))
dTransformed <- predict(preProcValues, d) %>% drop_na()

We drop the NA-rows (there are not too many), as further ahead these will be dropped anyway in bas.glm.

Note that not all participants participate in an equal number of experiments:

dTransformed %>% group_by(participant_ID) %>%  count() %>%
  ggplot(aes(x=participant_ID,y=n)) + geom_col()

We use Bayesian model averaging using the BAS-library to gain insight into which variables are important to be included. We force categorical variables to be fully included or not (hence we don’t allow certain factor levels to be included, while other levels of the same factors not). We insist on a model that contains angular_momentum, for else we cannot determine the required threshold value from the inferred model.

We use a uniform prior over all possible models and use the default mixture of Zellner-g-prior on the coefficients.

IT <- 20000 # nr of MCMC iterations
set.seed(15)
bma0 = bas.glm(outcome ~ ., data=dTransformed, method="MCMC", 
               MCMC.iterations=IT, family=binomial(), 
               include.always = ~ angular_momentum,
               modelprior=uniform(), force.heredity=TRUE)
print(summary(bma0),digits=2)
##                   P(B != 0 | Y)  model 1  model 2  model 3  model 4  model 5
## Intercept                  1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID16           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID19           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID22           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID6            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID5            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID12           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID2            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID15           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID20           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID11           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID21           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID1            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID10           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID23           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID24           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID13           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID7            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID14           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID4            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID17           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID8            1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID18           1.00    1.000    1.000    1.000    1.000    1.000
## participant_ID3            1.00    1.000    1.000    1.000    1.000    1.000
## gendermale                 0.13    0.000    0.000    0.000    0.000    0.000
## age                        0.17    0.000    0.000    0.000    0.000    0.000
## weight                     0.13    0.000    0.000    0.000    0.000    0.000
## length                     0.15    0.000    0.000    0.000    0.000    0.000
## skill_performance          0.90    1.000    1.000    1.000    1.000    1.000
## skill_effort               0.30    0.000    0.000    0.000    1.000    0.000
## velocity                   1.00    1.000    1.000    1.000    1.000    1.000
## pull_ID                    1.00    1.000    1.000    1.000    1.000    1.000
## pull_directionCW           0.67    1.000    0.000    1.000    1.000    1.000
## angular_momentum           1.00    1.000    1.000    1.000    1.000    1.000
## reaction_time              0.20    0.000    0.000    0.000    0.000    0.000
## phi_zero                   0.55    1.000    1.000    0.000    1.000    0.000
## phi_dot_zero               0.20    0.000    0.000    0.000    0.000    0.000
## delta_zero                 0.29    0.000    0.000    1.000    0.000    0.000
## delta_dot_zero             0.95    1.000    1.000    1.000    1.000    1.000
## psi_zero                   0.18    0.000    0.000    0.000    0.000    0.000
## Q_zero                     0.26    0.000    0.000    0.000    0.000    0.000
## BF                           NA    1.000    0.748    0.368    0.324    0.492
## PostProbs                    NA    0.031    0.027    0.023    0.017    0.016
## R2                           NA       NA       NA       NA       NA       NA
## dim                          NA   31.000   30.000   31.000   32.000   30.000
## logmarg                      NA -317.066 -317.356 -318.066 -318.192 -317.776
diagnostics(bma0) # should be on straight line

Experiments containing missing values are dropped in the analysis. As seen in the exploratory data analysis, the fraction of missing data is modest. The following figure visualises the models which get highest posterior probability. Black blocks denote that a variable is not contained in the model

image(bma0, rotate=F)

One way to choose a model, is by taking the model for which the marginal posterior inclusion probabilities exceed 1/2, the so called posterior median model.

aa <- tibble(probne0=bma0$probne0, names=bma0$namesx) 
partic_p <- aa %>% filter(names=="participant_ID1") %>% select(probne0)

p_inclusion <- aa %>% 
  filter(!str_detect(names,"participant_ID")) %>%
  filter(!str_detect(names,"Interc")) %>%
  add_row(names="participant_ID", probne0=as.numeric(partic_p)) %>% 
  mutate(names2= reorder(names, probne0, median))  %>% 
  ggplot(aes(x=probne0, y=names2)) + 
  geom_point() + geom_vline(xintercept=0.5, colour="red") +
  labs(x="posterior inclusion probability", y="")

p_inclusion

This suggests a model containing: angular_momentum, participant_ID, delta_ dot_ zero, pull_ID, velocity, skill_performance, pull_direction and phi_zero.

Now suppose we drop participant_ID from the list of predictors.

set.seed(13)
bma1 = bas.glm(outcome ~ ., data=dTransformed[,-1], method="MCMC", 
               include.always = ~ angular_momentum,
               MCMC.iterations=IT, family=binomial(), 
               modelprior=uniform(), force.heredity=TRUE)
print(summary(bma1),digits=2)
##                   P(B != 0 | Y)  model 1  model 2 model 3  model 4  model 5
## Intercept                  1.00    1.000    1.000    1.00    1.000    1.000
## gendermale                 0.99    1.000    1.000    1.00    1.000    1.000
## age                        1.00    1.000    1.000    1.00    1.000    1.000
## weight                     0.40    0.000    0.000    0.00    0.000    0.000
## length                     0.56    1.000    1.000    1.00    1.000    0.000
## skill_performance          0.40    0.000    1.000    0.00    1.000    0.000
## skill_effort               0.20    0.000    0.000    0.00    0.000    0.000
## velocity                   1.00    1.000    1.000    1.00    1.000    1.000
## pull_ID                    0.98    1.000    1.000    1.00    1.000    1.000
## pull_directionCW           0.54    0.000    0.000    1.00    1.000    0.000
## angular_momentum           1.00    1.000    1.000    1.00    1.000    1.000
## reaction_time              0.13    0.000    0.000    0.00    0.000    0.000
## phi_zero                   0.39    0.000    0.000    1.00    0.000    0.000
## phi_dot_zero               0.57    0.000    0.000    1.00    1.000    0.000
## delta_zero                 0.91    1.000    1.000    1.00    1.000    1.000
## delta_dot_zero             0.91    1.000    1.000    1.00    1.000    1.000
## psi_zero                   0.12    0.000    0.000    0.00    0.000    0.000
## Q_zero                     0.93    1.000    1.000    1.00    1.000    1.000
## BF                           NA    0.851    0.875    1.00    0.887    0.830
## PostProbs                    NA    0.027    0.021    0.02    0.019    0.018
## R2                           NA       NA       NA      NA       NA       NA
## dim                          NA   10.000   11.000   13.00   13.000    9.000
## logmarg                      NA -419.657 -419.629 -419.49 -419.615 -419.682
diagnostics(bma1) # should be on straight line 

image(bma1, rotate=F)

aa1 <- tibble(probne0=bma1$probne0, names=bma1$namesx) 
partic_p <- aa1 %>% filter(names=="participant_ID1") %>% select(probne0)

p_inclusion1 <- aa1 %>% 
  filter(!str_detect(names,"participant_ID")) %>%
  filter(!str_detect(names,"Interc")) %>%
  add_row(names="participant_ID", probne0=as.numeric(partic_p)) %>% 
  mutate(names2= reorder(names, probne0, median))  %>% 
  ggplot(aes(x=probne0, y=names2)) + 
  geom_point() + geom_vline(xintercept=0.5, colour="red") +
  labs(x="posterior inclusion probability", y="")

p_inclusion1

Not surprisingly, in this case some participant characteristics enter the median probability model.

We make a plot of confidence intervals for the coefficients, excluding those for participants.

c0 <- confint(coef(bma0))
c0_df <- tibble(coefficient=names(c0[,1]), mean=as.vector(c0[,3]),
                low=as.vector(c0[,1]), up= as.vector(c0[,2])) %>%
                mutate(type=str_detect(coefficient,"partic"))

c1 <- confint(coef(bma1))
c1_df <- tibble(coefficient=names(c1[,1]), mean=as.vector(c1[,3]), 
                low=as.vector(c1[,1]), up= as.vector(c1[,2])) %>%
                mutate(type=str_detect(coefficient,"partic"))

c01_df <- rbind(c0_df, c1_df) %>% filter(type==FALSE) %>%  
  mutate(participant_included=rep(c("yes", "no"),each=18))

c01_df %>% 
  ggplot(aes(x=coefficient, y=mean,color=participant_included)) +
  geom_point() + 
  geom_errorbar(aes(ymin=low, ymax=up)) +
  geom_hline(yintercept = 0)+ coord_flip() + 
  theme(legend.position = "bottom")

Multilevel approach

The preceding analysis has been helpful in identifying important predictors. It does however neglect the multilevel structure. It does not seem easy to specify that in bas.glm. For that reason, we refit the model, including the multilevel structure, using the brms package.

We first set the prior distributions on the coefficients to be zero-mean Normally distributed with standard-deviation of 2.

prior1 <- set_prior("normal(0, 2)", class = "b")

We fit a multilevel version based on the median probability model appearing from bma0.

set.seed(16)
fit0 <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
            skill_performance + delta_dot_zero + velocity +
            pull_ID + angular_momentum +  (1 | participant_ID), 
            data=dTransformed,
            family = bernoulli(link="logit"),
            prior=prior1,
            warmup = 5000, 
            iter = IT, 
            chains = 4, 
            seed = 123)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000167 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.67 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 1: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 1: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 1: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 1: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 1: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 8.15245 seconds (Warm-up)
## Chain 1:                31.6577 seconds (Sampling)
## Chain 1:                39.8101 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 7.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.76 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 2: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 2: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 2: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 2: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 2: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 7.94158 seconds (Warm-up)
## Chain 2:                32.8181 seconds (Sampling)
## Chain 2:                40.7597 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 6.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.65 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 3: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 3: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 3: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 3: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 3: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 7.77018 seconds (Warm-up)
## Chain 3:                30.4534 seconds (Sampling)
## Chain 3:                38.2236 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '1919117ad03e75562c036424e5d8ac5f' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 9e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.9 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 4: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 4: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 4: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 4: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 4: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 8.56809 seconds (Warm-up)
## Chain 4:                30.9096 seconds (Sampling)
## Chain 4:                39.4777 seconds (Total)
## Chain 4:
fit0 %>% summary()
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID) 
##    Data: dTransformed (Number of observations: 874) 
##   Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Group-Level Effects: 
## ~participant_ID (Number of levels: 24) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.69      0.62     2.67     5.07 1.00    11133    20445
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.48      0.72    -0.93     1.94 1.00     9167    15700
## phi_zero              0.42      0.18     0.07     0.77 1.00    41696    41662
## pull_directionCW     -0.49      0.25    -0.98    -0.01 1.00    47048    42765
## skill_performance     0.71      0.22     0.28     1.14 1.00    31549    39407
## delta_dot_zero       -0.49      0.15    -0.79    -0.18 1.00    43822    42088
## velocity              1.62      0.18     1.27     1.99 1.00    33574    39211
## pull_ID              -0.77      0.15    -1.06    -0.48 1.00    43260    40636
## angular_momentum      5.26      0.38     4.53     6.03 1.00    33461    36826
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
prior_summary(fit0)
##                 prior class              coef          group resp dpar nlpar
##          normal(0, 2)     b                                                 
##          normal(0, 2)     b  angular_momentum                               
##          normal(0, 2)     b    delta_dot_zero                               
##          normal(0, 2)     b         Intercept                               
##          normal(0, 2)     b          phi_zero                               
##          normal(0, 2)     b  pull_directionCW                               
##          normal(0, 2)     b           pull_ID                               
##          normal(0, 2)     b skill_performance                               
##          normal(0, 2)     b          velocity                               
##  student_t(3, 0, 2.5)    sd                                                 
##  student_t(3, 0, 2.5)    sd                   participant_ID                
##  student_t(3, 0, 2.5)    sd         Intercept participant_ID                
##  bound       source
##                user
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##             default
##        (vectorized)
##        (vectorized)

These coefficients are not exactly the same as under bma0 since

Compare Bayesian with frequentist results for model 0

dTransformed_reformat <- dTransformed %>%
                  mutate(outcomenum = as.numeric(outcome)-1)
glmer_fit0 <- glmer(outcomenum ~  phi_zero + pull_direction +
                      skill_performance + delta_dot_zero + velocity +
                      pull_ID + angular_momentum +  (1 | participant_ID), 
                     data=dTransformed_reformat, family = binomial)
glmer_fit0 %>% summary()
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: outcomenum ~ phi_zero + pull_direction + skill_performance +  
##     delta_dot_zero + velocity + pull_ID + angular_momentum +  
##     (1 | participant_ID)
##    Data: dTransformed_reformat
## 
##      AIC      BIC   logLik deviance df.resid 
##    624.1    667.0   -303.0    606.1      865 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6102 -0.2992 -0.0160  0.2774  7.3427 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  participant_ID (Intercept) 13.28    3.644   
## Number of obs: 874, groups:  participant_ID, 24
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.5533     0.7689   0.720  0.47180    
## phi_zero            0.4236     0.1795   2.360  0.01828 *  
## pull_directionCW   -0.5070     0.2500  -2.028  0.04257 *  
## skill_performance   0.7343     0.2252   3.261  0.00111 ** 
## delta_dot_zero     -0.4956     0.1561  -3.174  0.00150 ** 
## velocity            1.6616     0.1879   8.844  < 2e-16 ***
## pull_ID            -0.7837     0.1527  -5.132 2.87e-07 ***
## angular_momentum    5.3777     0.4015  13.394  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) phi_zr pll_CW skll_p dlt_d_ velcty pll_ID
## phi_zero    -0.033                                          
## pll_drctnCW -0.177  0.078                                   
## skll_prfrmn -0.023  0.049 -0.030                            
## delta_dt_zr  0.059 -0.258 -0.388 -0.014                     
## velocity     0.032  0.084 -0.098  0.477 -0.060              
## pull_ID      0.010 -0.033  0.023  0.267  0.125 -0.024       
## anglr_mmntm  0.049  0.120 -0.099  0.240 -0.205  0.608 -0.397

To facilitate comparison, we show both the Bayesian and frequentist fitted models:

sum0 <- fit0 %>% summary() %>% unclass() 
sum0$fixed
##                     Estimate Est.Error    l-95% CI    u-95% CI     Rhat
## Intercept          0.4829009 0.7246277 -0.93089505  1.93835361 1.000372
## phi_zero           0.4169201 0.1770367  0.07206348  0.76711668 1.000102
## pull_directionCW  -0.4937041 0.2457098 -0.97539265 -0.01218042 1.000055
## skill_performance  0.7059425 0.2203738  0.27837368  1.14201295 1.000031
## delta_dot_zero    -0.4865379 0.1541700 -0.78841881 -0.18386285 1.000090
## velocity           1.6197151 0.1828174  1.27426927  1.99044089 1.000170
## pull_ID           -0.7658691 0.1495807 -1.06476398 -0.47598085 1.000059
## angular_momentum   5.2557019 0.3801358  4.53472387  6.02727389 1.000292
##                    Bulk_ESS Tail_ESS
## Intercept          9166.944 15699.62
## phi_zero          41696.270 41662.04
## pull_directionCW  47047.902 42765.20
## skill_performance 31549.339 39407.11
## delta_dot_zero    43821.940 42087.76
## velocity          33574.146 39211.02
## pull_ID           43259.604 40636.39
## angular_momentum  33460.657 36825.72
sum0$random
## $participant_ID
##               Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 3.694535 0.6165921 2.673694  5.07287 1.0003 11133.41 20445.24

to be compared to

glmer0 <- glmer_fit0 %>% summary() %>% unclass() 
coefficients(glmer0)
##                     Estimate Std. Error    z value     Pr(>|z|)
## (Intercept)        0.5532895  0.7689346  0.7195534 4.718000e-01
## phi_zero           0.4236194  0.1795149  2.3598010 1.828474e-02
## pull_directionCW  -0.5070225  0.2500183 -2.0279420 4.256617e-02
## skill_performance  0.7342650  0.2251922  3.2606145 1.111711e-03
## delta_dot_zero    -0.4955544  0.1561312 -3.1739608 1.503740e-03
## velocity           1.6616009  0.1878872  8.8436099 9.267472e-19
## pull_ID           -0.7837200  0.1527248 -5.1315842 2.873137e-07
## angular_momentum   5.3777007  0.4015094 13.3937119 6.580788e-41
glmer0$varcor
##  Groups         Name        Std.Dev.
##  participant_ID (Intercept) 3.644

The posterior means from brm and estimates from glmer are very similar. Additionally, we obtain significance of all variables in the model.

Posterior predictive check

Posterior predictive checks are used to establish visually if the fitted model captures certain features in the data appropriately. Here we define such a check.

The function checks2 compares the fraction of falls for the total number of perturbations of each participant to that based on 1000 predictions of the fitted model. It does so separately for each of the velocities, which are in the set {6,12,18}

checks2 <- function(d, fit) # fraction falling by participant and factor(velocity)
{
  pred <- posterior_predict(fit, ndraws=1000)
  fracfalling1 <- d %>% 
    mutate(velocity=as.factor(round(velocity,1))) %>%     
    mutate(predicted_average=colMeans(pred)) %>% 
    group_by(velocity, participant_ID) %>%
    summarise(n=n(), observed=mean(outcome=="1"),
              predicted=mean(predicted_average)) %>%
    pivot_longer(cols=c(predicted, observed), names_to="type", values_to="fraction")
  
  p <- fracfalling1 %>% 
    ggplot() + 
    geom_point(aes(x=velocity, y = fraction,colour=type)) +
    labs(x="velocity (standardized)", y="fraction of falls of the total number of perturbations") 
  p + facet_wrap(~participant_ID) 
}

checks2(dTransformed, fit0)

Prior sensitivity

We check the effect of the specified prior distribution on the estimates for the coefficients.

prior_vague <- set_prior("normal(0, 5)", class = "b") # way larger sd

prior_centred  <- set_prior("normal(0, 1)", class = "b") 
# somewhat small sd, possibly pulls estimates towards zero.

Also fit the model with these priors:

fit0_vague <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
            skill_performance + delta_dot_zero + velocity +
            pull_ID + angular_momentum +  (1 | participant_ID), 
            data=dTransformed,
            family = bernoulli(link="logit"),
            prior=prior_vague,
            warmup = 5000, 
            iter = IT, 
            chains = 4, 
            seed = 123)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000224 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.24 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 1: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 1: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 1: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 1: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 1: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 8.72305 seconds (Warm-up)
## Chain 1:                32.1648 seconds (Sampling)
## Chain 1:                40.8878 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 8.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.83 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 2: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 2: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 2: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 2: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 2: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 8.37378 seconds (Warm-up)
## Chain 2:                32.1127 seconds (Sampling)
## Chain 2:                40.4865 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 5.8e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 3: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 3: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 3: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 3: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 3: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 9.01871 seconds (Warm-up)
## Chain 3:                33.0101 seconds (Sampling)
## Chain 3:                42.0288 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'ce805167e8faadb918648a9d001dabfe' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.6 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 4: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 4: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 4: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 4: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 4: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 8.67467 seconds (Warm-up)
## Chain 4:                31.6836 seconds (Sampling)
## Chain 4:                40.3583 seconds (Total)
## Chain 4:
fit0_centred <- brm(outcome ~ 0 + Intercept + phi_zero + pull_direction +
            skill_performance + delta_dot_zero + velocity +
            pull_ID + angular_momentum +  (1 | participant_ID), 
            data=dTransformed,
            family = bernoulli(link="logit"),
            prior=prior_centred,
            warmup = 5000, 
            iter = IT, 
            chains = 4, 
            seed = 123)
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## 
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0.000155 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.55 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 1: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 1: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 1: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 1: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 1: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 1: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 1: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 1: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 1: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 1: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 1: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 7.19273 seconds (Warm-up)
## Chain 1:                26.2807 seconds (Sampling)
## Chain 1:                33.4734 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 6.9e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 2: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 2: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 2: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 2: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 2: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 2: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 2: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 2: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 2: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 2: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 2: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 6.8328 seconds (Warm-up)
## Chain 2:                24.6561 seconds (Sampling)
## Chain 2:                31.4889 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 5.9e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 3: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 3: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 3: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 3: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 3: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 3: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 3: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 3: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 3: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 3: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 3: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 7.05009 seconds (Warm-up)
## Chain 3:                23.7962 seconds (Sampling)
## Chain 3:                30.8463 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'a07710f062e3c373f1c44960e5a29b2c' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 6.1e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.61 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:     1 / 20000 [  0%]  (Warmup)
## Chain 4: Iteration:  2000 / 20000 [ 10%]  (Warmup)
## Chain 4: Iteration:  4000 / 20000 [ 20%]  (Warmup)
## Chain 4: Iteration:  5001 / 20000 [ 25%]  (Sampling)
## Chain 4: Iteration:  7000 / 20000 [ 35%]  (Sampling)
## Chain 4: Iteration:  9000 / 20000 [ 45%]  (Sampling)
## Chain 4: Iteration: 11000 / 20000 [ 55%]  (Sampling)
## Chain 4: Iteration: 13000 / 20000 [ 65%]  (Sampling)
## Chain 4: Iteration: 15000 / 20000 [ 75%]  (Sampling)
## Chain 4: Iteration: 17000 / 20000 [ 85%]  (Sampling)
## Chain 4: Iteration: 19000 / 20000 [ 95%]  (Sampling)
## Chain 4: Iteration: 20000 / 20000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 7.07604 seconds (Warm-up)
## Chain 4:                22.6028 seconds (Sampling)
## Chain 4:                29.6788 seconds (Total)
## Chain 4:
fit0 %>% summary()
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID) 
##    Data: dTransformed (Number of observations: 874) 
##   Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Group-Level Effects: 
## ~participant_ID (Number of levels: 24) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.69      0.62     2.67     5.07 1.00    11133    20445
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.48      0.72    -0.93     1.94 1.00     9167    15700
## phi_zero              0.42      0.18     0.07     0.77 1.00    41696    41662
## pull_directionCW     -0.49      0.25    -0.98    -0.01 1.00    47048    42765
## skill_performance     0.71      0.22     0.28     1.14 1.00    31549    39407
## delta_dot_zero       -0.49      0.15    -0.79    -0.18 1.00    43822    42088
## velocity              1.62      0.18     1.27     1.99 1.00    33574    39211
## pull_ID              -0.77      0.15    -1.06    -0.48 1.00    43260    40636
## angular_momentum      5.26      0.38     4.53     6.03 1.00    33461    36826
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
fit0_vague %>% summary()
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID) 
##    Data: dTransformed (Number of observations: 874) 
##   Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Group-Level Effects: 
## ~participant_ID (Number of levels: 24) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.84      0.65     2.78     5.30 1.00    11541    20586
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.57      0.80    -0.99     2.14 1.00     9069    15652
## phi_zero              0.43      0.18     0.08     0.79 1.00    41898    40680
## pull_directionCW     -0.51      0.25    -1.01    -0.03 1.00    46950    40973
## skill_performance     0.74      0.23     0.31     1.19 1.00    31771    37675
## delta_dot_zero       -0.50      0.16    -0.81    -0.20 1.00    44348    42126
## velocity              1.68      0.19     1.33     2.06 1.00    32655    36795
## pull_ID              -0.79      0.15    -1.10    -0.50 1.00    41567    39251
## angular_momentum      5.46      0.40     4.71     6.28 1.00    31567    34596
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
fit0_centred %>% summary()
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: outcome ~ 0 + Intercept + phi_zero + pull_direction + skill_performance + delta_dot_zero + velocity + pull_ID + angular_momentum + (1 | participant_ID) 
##    Data: dTransformed (Number of observations: 874) 
##   Draws: 4 chains, each with iter = 20000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Group-Level Effects: 
## ~participant_ID (Number of levels: 24) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     3.32      0.55     2.41     4.57 1.00    11275    22038
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept             0.33      0.58    -0.82     1.45 1.00     8446    15627
## phi_zero              0.38      0.17     0.05     0.71 1.00    36554    40778
## pull_directionCW     -0.44      0.23    -0.89     0.01 1.00    44925    42130
## skill_performance     0.61      0.20     0.22     1.02 1.00    29294    37257
## delta_dot_zero       -0.44      0.15    -0.73    -0.16 1.00    43493    44123
## velocity              1.45      0.17     1.13     1.78 1.00    32566    37339
## pull_ID              -0.69      0.14    -0.97    -0.42 1.00    41114    39827
## angular_momentum      4.73      0.32     4.13     5.38 1.00    31855    39929
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

From here, we see that the differences with the “vague” prior are small; the “centred” prior shrinks coefficients more towards zero and appears to have some effect on the estimates.

Visualise effect of numerical predictors

Some example plots to illustrate the effect of a particular variable on outcome.

conditional_effects(fit0, effects = "pull_direction")

conditional_effects(fit0, effects = c("velocity"))

conditional_effects(fit0, effects = c("delta_dot_zero"))

conditional_effects(fit0, effects = c("angular_momentum:velocity"))

Visualising the posterior

We visualise the posterior by showing density- and traceplots of posterior samples.

fit0 %>% plot(newpage=FALSE, ask=FALSE)

We extract IT posterior samples from the fitted model.

ps <- brms::as_draws_df(fit0) %>% tibble() %>% sample_n(IT)

We plot all parameters, except for participant effects

ln <- c("b_Intercept", "b_skill_performance", "b_velocity", "b_phi_zero",
        "b_pull_ID", "b_pull_directionCW",  "b_angular_momentum",
         "b_delta_dot_zero","sd_participant_ID__Intercept")

ps_coef <- ps %>% dplyr::select(-contains("r_particip")) %>%
  dplyr::select(contains("b_") | contains("sd_")) %>%
  tidyr::pivot_longer(cols=ln,names_to="coefficient", values_to="value") %>%
  mutate(coefficient=str_remove(coefficient,"b_")) %>% 
  mutate(coefficient = recode(coefficient, `sd_participant_ID__Intercept` = "sd_participant"))

ps_coef %>% ggplot(aes(x=value, y = reorder(coefficient, value,median))) +
    stat_halfeye( point_size=1, interval_size=0.5)+
    labs(x="coefficient", y = "")

Similarly for participant effects

ps_sub <- ps %>% dplyr::select(contains("r_particip"))
ps_sublong <-  ps_sub %>% tidyr::pivot_longer(cols=contains("r_particip"),
                names_to="participant", values_to="value") %>%
                mutate(participant_ID=as.factor(readr::parse_number(participant)))
ps_complete <- left_join(ps_sublong, d, by ="participant_ID")

ps_complete %>% ggplot(aes(x=value, 
            y = reorder(participant_ID, value,median),fill=age )) +
            stat_halfeye(point_colour="orange", point_size=1, interval_size=0.5)+
         #   labs(x="coefficient", y = "participant_ID")
    labs(x=TeX(r'($\beta_i$)'), y="participant_ID")

rm(ps_coef, ps_sub, ps_sublong, ps_complete) # free memory

Computing perturbation threshold for a few participants in the study.

Consider the results for following persons:

pIDs <- c(1,1,9,9,11,11,11,12,12) # arbitrarily choose a few participants
d1 <- d %>% filter(participant_ID==as.character(1)) %>% slice(1)
d9 <- d %>% filter(participant_ID==as.character(9)) %>% slice(1)
d11 <- d %>% filter(participant_ID==as.character(11)) %>% slice(1)
d12 <- d %>% filter(participant_ID==as.character(12)) %>% slice(1)
d_existing <- bind_rows(d1,d1,d9,d9,d11,d11,d11,d12,d12) %>% 
  mutate(participant_ID2=c("P1_v6","P1_v12","P9_v6", "P9_v12", "P11_v6","P11_v12","P11_v18", "P12_v6", "P12_v12" )) %>%
  mutate(velocity=c(6,12,6,12,6,12,18,6,12)) %>%
  mutate(phi_zero=rep(0,9), pull_direction=rep("CW",9), delta_dot_zero=rep(0,9), pull_ID=rep(1,9))

d_existingTransformed <- predict(preProcValues, d_existing)


nps <- nrow(ps) # nr of posterior samples
zz = rep(0,nps)
z <- data.frame("P1_v6" = zz,"P1_v12" = zz,"P9_v6" = zz, "P9_v12" = zz, "P11_v6" = zz,
                "P11_v12" = zz,"P11_v18" = zz, "P12_v6" = zz, "P12_v12" = zz)
head(z)
##   P1_v6 P1_v12 P9_v6 P9_v12 P11_v6 P11_v12 P11_v18 P12_v6 P12_v12
## 1     0      0     0      0      0       0       0      0       0
## 2     0      0     0      0      0       0       0      0       0
## 3     0      0     0      0      0       0       0      0       0
## 4     0      0     0      0      0       0       0      0       0
## 5     0      0     0      0      0       0       0      0       0
## 6     0      0     0      0      0       0       0      0       0

Function to compute threshold function

angular_threshold <- function(participant_ch, postdraw, participant_postdraw)
{
    xx <- participant_postdraw +
        postdraw$b_Intercept + 
        postdraw$b_skill_performance * participant_ch$skill_performance +
        postdraw$b_velocity * participant_ch$velocity+
        postdraw$b_pull_ID * participant_ch$pull_ID +
        postdraw$b_pull_directionCW * (participant_ch$pull_direction=="CW") + 
      #  postdraw$b_phi_zero * participant_ch$phi_zero +
        postdraw$b_delta_dot_zero * participant_ch$delta_dot_zero 
    out <- -xx/(postdraw$b_angular_momentum) 
    out
}  

for (j in 1:9)
{
  ID <- pIDs[j]
  # select for the participant posterior draws from its intercept:
  lab <- paste0("r_participant_ID[",as.character(ID),",Intercept]")
  raneff <- ps[,lab] %>% deframe()
  #  choose participant experimental setting:
  participant_ch <- d_existingTransformed[j,] 

  for (i in 1:nps)  # for each posterior sample
  { 
    postdraw <- ps %>% slice(i) 
    z[i, j] <- angular_threshold(participant_ch, postdraw, raneff[i])
  }
}
m_ang <- mean(na.omit(d)$angular_momentum)
sd_ang <- sd(na.omit(d)$angular_momentum)

angular_threshold_df <- apply(z,2, function(x) {  m_ang + sd_ang * x  }  )  # rescale each column
angular_threshold_tidy <- as_tibble(angular_threshold_df) %>%
                          pivot_longer(cols=1:9, names_to="participant", values_to="eta")

Plotting the results

angular_threshold_tidy %>% 
  #ggplot(aes(x=eta, y=reorder(participant, eta, median))) +
  ggplot(aes(x=eta, y=participant)) +
  stat_halfeye(fill="lightblue") +
  labs(x=TeX(r'($\Delta L_{50\%}$)'), y ="") +
  ggtitle("Threshold value for angular momentum")

Computing perturbation threshold for six new fictive participants.

We manually construct “new” participants, for which we will compute the perturbation threshold (for the variables in model 5 we provide values, while all remaining columns are just filled with NA-values).

New fictive participants

dpred <- tibble(participant_ID=c("new_v12_sp40", "new_v12_sp90", "new_v6_sp40", "new_v6_sp90","new_v18_sp40","new_v25_sp40"),
                velocity=c(12,12,6,6,18,25), age = rep(NA,6),
                skill_performance=c(40,90,40,90,40,40), pull_ID=rep(1,6),
                pull_direction=rep("CW",6), weight=rep(NA,6), length=rep(NA,6),
                skill_effort = rep(NA,6), initial_search=rep(NA,6), phi_zero=rep(0,6),
                phi_dot_zero=rep(NA,6), delta_zero=rep(NA,6), delta_dot_zero=rep(0,6),
                psi_zero=rep(NA,6), Q_zero=rep(NA,6), comments=rep(NA,6),
                exclude=rep(NA,6), reaction_time=rep(NA,6), pull_force=rep(NA,6),
                angular_momentum=rep(NA,6))

Now we will predict for these new participants:

dpred_tf <- predict(preProcValues, dpred)  # tf=transformed as we did for the model
dpred_tf$delta_dot_zero <- rep(0,6)

nparticipants = nrow(dpred)  # nr of participants for which we wish to predict
zz = rep(0,nps)
z <- data.frame(new_participant1 = zz,  new_participant2 = zz,
                new_participant3 = zz,new_participant4 = zz,
                new_participant5 = zz, new_participant6 = zz)

for (j in 1:nparticipants)  # for each participant we wish to find threshold value
{
  for (i in 1:nps)  # for each participant
  { 
    particip =  rnorm(1,sd=ps$sd_participant_ID__Intercept[i])
    z[i, j] <- angular_threshold(dpred_tf[j,], ps[i,], particip)
  }
}  
# rescale each column:
angular_threshold_df_new <- apply(z,2, function(x) { m_ang + sd_ang * x})  
angular_threshold_tidy_new <- as_tibble(angular_threshold_df_new) %>%
      pivot_longer(cols=contains("new"), names_to="participant", values_to="eta") %>%
      mutate(participant=fct_recode(participant, 
                                    s40v12="new_participant1",
                                    s90v12="new_participant2",
                                    s40v6="new_participant3",
                                    s90v6="new_participant4",
                                    s40v18="new_participant5",
                                    s40v25="new_participant6"))

Plotting the results

angular_threshold_tidy_new %>% 
#  ggplot(aes(x=eta, y=reorder(participant, eta, median))) +
  ggplot(aes(x=eta, y=participant)) +
  stat_halfeye(fill="lightblue") +
  labs(x=TeX(r'($\Delta L_{50\%}$)'), y ="") +
  ggtitle("Threshold value for angular momentum")

We can also make a combined plot for participants which were part of the study and “new” participants:

nr <- nrow(angular_threshold_tidy)
nr_new <- nrow(angular_threshold_tidy_new)
factor_new <- c(rep("no",nr), rep("yes",nr_new))

bind_rows(angular_threshold_tidy, angular_threshold_tidy_new) %>% 
    mutate(new=factor_new) %>%  
    ggplot(aes(x=eta, y=reorder(participant, eta, median), fill=new)) +
    stat_halfeye(fill="lightblue") + labs(x="eta", y = "participant")

This shows that we can fairly accurately find the threshold value for participants within the study, but have large uncertainty about this value for “future” unseen participants. This is somewhat natural, as we don’t know the inherent cycling performance of unseen participants.